Test 3: Caso 2D - Cavidad

En el flujo con una cavidad rectangular, el fluido se encuentra confinado en una región rectangular de dimensiones $L\times W$ donde el movimiento se debe al desplazamiento transversal de una de las caras. En este sentido, la velocidad de la cara superior es de $U_0$ y el de las demás interfaces es de $0$. De acuerdo con esto, se tienen las siguientes consideraciones:

  • El flujo se presenta únicamente en las direcciones $x$ y $y$ ($w=0$).
  • El flujo no depende de la coordenada $z$ ($\frac{\partial}{\partial z}=0$).
  • El fluido se encuentra en estado estacionario ($\frac{\partial}{\partial t}=0$).
  • No hay fuerzas externas afectando el flujo ($f_x=f_y=f_z=0$).

De aquí, las ecuaciones de continuidad y de Navier-Stokes se reducen a:

\begin{align} \frac{\partial \rho}{\partial t} + \nabla\cdot(\rho \vec{u})=0 \quad &\to\quad \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0 \\ \rho\left(\frac{\partial \vec{u}}{\partial t}+\vec{u}\cdot\nabla\vec{u}\right)=\rho\vec{f}-\nabla p+\mu\nabla^2\vec{u}\quad &\to\quad \begin{cases} \rho \left(u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)=-\frac{\partial p}{\partial x}+\mu\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}\right) \\ \rho \left(u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}\right)=-\frac{\partial p}{\partial v}+\mu\left(\frac{\partial^2v}{\partial x^2}+\frac{\partial^2v}{\partial y^2}\right) \end{cases} \end{align}

Se puede ver entonces que la ecuación resultante es un caso particular del problema de convección-difusión estable:

$$\frac{\partial f}{\partial t} + \vec{u}\cdot\nabla f=k\nabla^2f+\frac{s}{\rho c_p}$$

Donde $\frac{\partial f}{\partial t}=0$, $k = \frac{\mu}{\rho}=\nu$, y $\frac{s}{c_p}=-\frac{\partial p}{\partial *}$ para $f=u, *=x$ y $f=v, *=y$, y las condiciones de frontera corresponden a valores conocidos (tipo Dirichlet).

Ahora, para aplicar el método de elementos espectrales, es necesario realizar los siguientes pasos:

  • Introducción de la información del problema
  • Definición de los elementos a usar (cantidad $N_E$ y distribución, así como sus órdenes $N_P(i)$)
    • Selección de los nodos de interpolación a usar
    • Generación de la matriz de conectividad
  • Generación de la matriz del problema
    • Generación de la matriz de difusión
      • Calculo de las funciones de interpolación
      • Calculo de las derivadas de las funciones de interpolación
      • Calculo de la función de transformación y sus derivadas
      • Calculo de los gradientes de las funciones de interpolación
      • Calculo del factor de corrección (cambio de base) $h_s$
    • Generación de las matrices de presión
    • Generación del vector b
      • Implementación de las condiciones de frontera
  • Solución del sistema resultante
  • Gráfica de la solución obtenida

Cada uno de los pasos se tratará en detalle a continuación


In [1]:
# Prepara el ambiente de trabajo con las librerías a utilizar

import numpy as np               # Ofrece funciones para operaciones básicas con arreglos
import scipy.linalg as lin       # Ofrece funciones para operaciones de álgebra lineal
import matplotlib.pyplot as plt  # Permite incluir un ambiente de visualización
from mpl_toolkits.mplot3d.axes3d import Axes3D
import timeit                    # Permite la medición de tiempos para los distintos algoritmos
from sem2D import *               # Agrupa las funciones externas en un archivo
# from __future__ import division  # Corrige posibles problemas que surgen a la hora de dividir enteros
#%pylab inline

Introducción de la información del problema

Se introducen los parámetros del problema. En este caso se tienen los parámetros $L$, $W$ y $U_0$.


In [2]:
## Valores de prueba

L = 2.0
W = 2.0
mu = 1.
U0 = 1.

Definición de los elementos a usar

Se define primero el número de elementos y su distribución, así como los órdenes de cada uno de los mismos. Por simplicidad, inicialmente se toman elementos uniformemente espaciados.


In [3]:
Ne = np.array([12,8])   # Indica la cantidad de elementos a lo largo de los ejes x y y
Np = np.array([3,2])   # Indica la cantidad de divisiones sobre cada elemento en xi y en eta

# Genera los puntos correspondientes a los vertices de cada elemento
(px, hx) = np.linspace(-L/2,L/2,Ne[0]+1,retstep=True)
(py, hy) = np.linspace(-W/2,W/2,Ne[1]+1,retstep=True)

Selección de los nodos de interpolación a usar

Conociendo el número de elementos y su distribución, es posible generar los nodos correspondientes a los vértices de cada elemento. Sin embargo, como no se trabaja con elementos lineales, es necesario generar los nodos de interpolación para cada elemento. Para ello, es necesario definir la familia de polinomios sobre la cual van a estar basados estos nodos. En este caso se utilizan los nodos de interpolación de Lobatto.

def genNodes(Ne,Np,px,py,ordn='xy'): #================================================= # Genera los nodos de interpolacion segun el # numero de elementos y sus ordenes. Esta es la # numeracion local de los nodos: (l,i) # # Ne = Np.size # max(Np) < 6 #================================================= ne = Ne[0]*Ne[1] # Recupera el numero de elementos npe = (Np[0]+1)*(Np[1]+1) # Recupera el numero de nodos por elemento x = np.zeros([ne,npe]) # Coordenadas x por elemento y = np.zeros([ne,npe]) # Coordenadas y por elemento # Genera los nodos de interpolacion de Lobatto dx = qlobatto(Np[0]-1,'z') dy = qlobatto(Np[1]-1,'z') for l in range(ne): cx = l % Ne[0] cy = l / Ne[0] midx = 0.5 * (px[cx+1]+px[cx]) disx = 0.5 * (px[cx+1]-px[cx]) midy = 0.5 * (py[cy+1]+py[cy]) disy = 0.5 * (py[cy+1]-py[cy]) pxe = midx+disx*dx pye = midy+disy*dy (xv,yv) = np.meshgrid(pxe,pye,indexing=ordn) x[l,:] = xv.reshape(npe) y[l,:] = yv.reshape(npe) return x, y, ne, npe

In [4]:
(x, y, ne, npe) = genNodes(Ne,Np,px,py)

# Grafica la malla de puntos
for l in range(ne):
    plt.plot([x[l,0],x[l,Np[0]],x[l,-1],x[l,-Np[0]-1],x[l,0]],
             [y[l,0],y[l,Np[0]],y[l,-1],y[l,-Np[0]-1],y[l,0]],
             'r--')
    plt.hold(True)
plt.plot(x,y,'bo',ms=4)
plt.xlabel('coordenada x')
plt.ylabel('coordenada y')
plt.show()

Generación de la matriz de conectividad

Una vez se definen los nodos, se utiliza la matriz de conectividad para almacenar de manera ordenada la numeración local y global de los mismos

def genConn(Ne,Np,ne,npe,x,y,h): #================================================= # Genera los vectores de numeracion global, la # matriz de conectividad y el vector de frontera # # max(Np) < 6 #================================================= # Calcula el numero de nodos globales para definir el tamano de las matrices contenedoras ng = ne*npe - ((Ne[0]-1)*(Np[1]+1)+(Ne[1]-1)*(Np[0]+1)+(Ne[0]-1)*(Ne[1]-1)*(sum(Np)+1)) coords = np.zeros([ng,2]) # Almacena las coordenadas de todos los nodos en un vector sin redundancia C = np.zeros((ne,npe),int) # Matriz de conectividad gfl = np.zeros([ng,3]) # Indicador de nodos de frontera [indicador, vel. x, vel. y] # Se introduce manualmente el primer nodo coords[0,0] = x[0,0] coords[0,1] = y[0,0] # C[0,0] = 0 Esta condicion ya se tiene de la inicializacion de C gfl[0,0] = 1 # Indica un nodo de frontera (el valor de la velocidad es de [0,0]) cont = 1 # Contador global tol = h*1e-6 # Tolerancia para identificacion de nodos iguales for l in range(ne): for k in range(npe): existe = False; for n in range(cont): if(np.abs(coords[n,0]-x[l,k])+np.abs(coords[n,1]-y[l,k]) < tol): existe = True C[l,k] = n # El nodo (l,m) corresponde al nodo global n if(not existe): # Registra el nodo si no existe todavía coords[cont,0] = x[l,k] coords[cont,1] = y[l,k] C[l,k] = cont # Identifica los nodos de frontera y actualiza el indicador if (L/2.-np.abs(x[l,k]) < tol or W/2.-np.abs(y[l,k]) < tol): gfl[cont,0] = 1 # Si se encuentra en la frontera superior actualiza la velocidad x if(W/2.-y[l,k] < tol): gfl[cont,1] = U0 cont += 1 return coords, C, gfl, ng

In [5]:
h = min(hx,hy)
(coords, C, gfl, ng) = genConn(Ne,Np,ne,npe,x,y,L,W,U0,h,'quad')

Generación de la matriz del problema

Para la generación de la matriz del problema es necesario definir las matrices de difusión, masa y advección por cada elemento:

$$A_{ij}^{(\ell)}=\iint_{E_\ell} \nabla\psi_i^{(\ell)}\cdot\nabla\psi_j^{(\ell)}dxdy \quad B_{ij}^{(\ell)}=\iint_{E_\ell} \psi_i^{(\ell)}\psi_j^{(\ell)}dxdy \quad C_{ij}^{(\ell)}=\iint_{E_\ell} \psi_i^{(\ell)}\vec{u}\cdot\nabla\psi_j^{(\ell)}dxdy$$

Estas integrales se obtienen utilizando la cuadratura de Lobatto sobre cada una de las dimensiones, ya que las funciones de interpolación son bilineales y se escriben de la forma:

$$ \psi_i(\xi,\eta) = \mathcal{L}_i(\xi)\mathcal{W}_i(\eta) $$

donde $\mathcal{L}_i(\xi)$ y $\mathcal{W}_i(\eta)$ son las funciones de interpolación (unidimensionales) correspondientes. Por simplicidad, se trabaja sobre el elemento estandar con coordenadas $(-1,-1),(1,-1),(-1,1) \text{ y }(1,1)$ y se alude a un cambio de coordenadas.

Generación de la matriz de difusión

Comenzamos con la matriz que define el sistema de ecuaciones, o la matriz de difusión. Dado que esta depende de los gradientes de las funciones de interpolación, es necesario realizar todos los pasos de construcción de los términos como sigue.

Calculo de las funciones de interpolación

En este caso, debido al orden variable del método aplicado, se generan las funciones de interpolación a partir de los polinomios de interpolación de Lagrange. Para simplificar la notación se considera la biyección: $\{0,1,\dots,Np\} \to \{0,\dots,m\}\times\{0,\dots,n\}$ definida por $i \mapsto (j,k)$

$$ \psi_i(\xi,\eta) = \psi_{(j,k)}(\xi,\eta) = \mathcal{L}_j(\xi)\mathcal{W}_k(\eta) = \left(\prod_{l\neq j}^m\frac{\xi-\xi_l}{\xi_j-\xi_l}\right) \left(\prod_{l\neq k}^n\frac{\eta-\eta_l}{\eta_k-\eta_l}\right)$$

Calculo de las derivadas de las funciones de interpolación

De acuerdo con esta descripción, se tiene entonces que las derivadas con respecto a $\xi$ y $\eta$ son entonces:

\begin{align} \frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\xi} &= \psi_{(j,k)}(\xi,\eta) \sum_{l\neq j}^m \frac{1}{\xi-\xi_l} = \left[\sum_{l\neq j}^m \left(\frac{1}{\xi_j-\xi_l}\prod_{p\neq j,l}^m \frac{\xi-\xi_p}{\xi_j-\xi_p}\right)\right]\left(\prod_{l\neq k}^n\frac{\eta-\eta_l}{\eta_k-\eta_l}\right) \\ \frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\eta} &= \psi_{(j,k)}(\xi,\eta) \sum_{l\neq k}^m \frac{1}{\eta-\eta_l} = \left(\prod_{l\neq j}^n\frac{\xi-\xi_l}{\xi_j-\xi_l}\right) \left[\sum_{l\neq k}^m \left(\frac{1}{\eta_k-\eta_l}\prod_{p\neq k,l}^m \frac{\eta-\eta_p}{\eta_k-\eta_p}\right)\right] \end{align}

Calculo de la función de transformación y sus derivadas

De este modo, la función de transformación de coordenadas corresponde a:

$$ \vec{x} = \sum_{i=0}^{N_p-1} \vec{x}_i^E \psi_i(\xi,\eta) = \sum_{i=0}^{N_p-1} \vec{x}_i^E \psi_{(j,k)}(\xi,\eta) = \sum_{i=0}^{N_p-1} \vec{x}_i^E \mathcal{L}_j(\xi)\mathcal{W}_k(\eta) $$

y por lo tanto sus derivadas son:

\begin{align} \frac{\partial x}{\partial \xi} &= \sum_{i=0}^{N_p-1}x_i^E\frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\xi} &\quad \frac{\partial x}{\partial \eta} &= \sum_{i=0}^{N_p-1}x_i^E\frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\eta} \\ \frac{\partial y}{\partial \xi} &= \sum_{i=0}^{N_p-1}y_i^E\frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\xi} &\quad \frac{\partial y}{\partial \eta} &= \sum_{i=0}^{N_p-1}y_i^E\frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\eta} \end{align}

Calculo de los gradientes de las funciones de interpolación

Ya con estos valores es posible soluciones los gradientes de las funciones solucionando los sistemas:

$$ \begin{pmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{pmatrix} \nabla\psi_i = {\frac{\partial\psi_i}{\partial\xi} \choose \frac{\partial\psi_i}{\partial\xi}} $$

Calculo del factor de corrección (cambio de base) $h_s$

De manera que el coeficiente de cambio de coordenadas corresponde a:

$$ Det[J] = \begin{vmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{vmatrix} = \frac{\partial x}{\partial \xi}\frac{\partial y}{\partial \eta}-\frac{\partial x}{\partial \eta}\frac{\partial y}{\partial \xi} $$

Generación de las matrices de presión

Ahora bien, como tampoco se conocen los valores de presión, es necesario incorporarlas también como incógnitas. Para ello se calculan los vectores: $$ D_{i}^{(\ell)}=\iint_{E_\ell} \nabla\psi_i^{(\ell)}dxdy$$

Generación del vector b

Para el otro lado de la ecuación se implementa la función de entrada. Para esto se recorre sobre todos los elementos y se van calculando los valores para la función de manera correspondiente. Estos valores se acomodan al sistema global aprovechando la matriz de conectividad.

Implementación de las condiciones de frontera

Así mismo, es necesario modificar la matriz global de manera que se incorporen las condiciones de frontera de Dirichlet. Si $l$ es un nodo de frontera, esto se hace a través de los pasos:

  • Sustituir las entradas $b_i$ por $b_i-D_{il}f_l$
  • Cambiar el término $b_l$ por $f_l$
  • Igualar las entradas de la fila $l$ y columna $l$ por 0 en la matriz global
  • Reajustar el término $D_{ll}$ a $1$.
def matDiff(Ne,Np,ne,npe,ng,x,y,C,gfl,mu=1): #================================================= # Genera la matriz de difusion/rigidez y el vector # del lado derecho del sistema sujeto a condicio- # nes de Dirichlet # # max(Np) < 6 #================================================= # Recupera los nodos de interpolacion de Lobatto dx = qlobatto(Np[0]-1,'z')zpres = solp-np.mean(solp) dy = qlobatto(Np[1]-1,'z')zpres = solp-np.mean(solp) (xi,wx) = qlobatto(Np[0],'b') # Utiliza un orden mayor para tener una integral exacta (eta,wy) = qlobatto(Np[1],'b') # Utiliza un orden mayor para tener una integral exacta NQ = (Np[0]+2)*(Np[1]+2) # Numero de nodos a usar en la cuadratura (ix,iy) = np.meshgrid(range(Np[0]+2),range(Np[1]+2),indexing='xy') M_ind = np.array([ix.reshape(NQ),iy.reshape(NQ)]).T # Matriz de indices para la cuadratura Mat_dif = np.zeros([2*ng+ne,2*ng+ne]) # Inicializa la matriz de difusion global for l in range(ne): elm_dm = np.zeros([npe,npe]) # Inicializa la matriz de difusion elemental vecD = np.zeros([npe,2]) # Realiza la cuadratura de Gauss for nq in M_ind: psi = np.ones(npe) # Almacena las funciones de interpolacion dpsixi = np.zeros(npe) # Almacena las derivadas con respecto a xi dpsieta = np.zeros(npe) # Almacena las derivadas con respecto a eta grpsi = np.zeros([npe,2]) # Almacena los gradientes de las funciones de interpolacion dxdxi = 0. dxdeta = 0. dydxi = 0. dydeta = 0. for i in range(npe): #Realiza los calculos sobre cada elemento divx = i % (Np[0]+1) divy = i / (Np[0]+1) # Calcula las funciones de interpolacion fx = 1 for ind in range(Np[0]+1): if(ind != divx): fx *= (xi[nq[0]]-dx[ind])/(dx[divx]-dx[ind]) # print fx, fy = 1 for ind in range(Np[1]+1): if(ind != divy): fy *= (eta[nq[1]]-dy[ind])/(dy[divy]-dy[ind]) psi[i] = fx*fy # print fy, psi[i] # Calcula las derivadas con respecto a las dos variables for ind in range(Np[0]+1): if(ind != divx): prod = 1 for ind2 in range(Np[0]+1): if(ind2 != divx and ind2 != ind): prod *= (xi[nq[0]]-dx[ind2])/(dx[divx]-dx[ind2]) dpsixi[i] += 1./(dx[divx]-dx[ind])*prod dpsixi[i] *= fy for ind in range(Np[1]+1): if(ind != divy): prod = 1 for ind2 in range(Np[1]+1): if(ind2 != divy and ind2 != ind): prod *= (eta[nq[1]]-dy[ind2])/(dy[divy]-dy[ind2]) dpsieta[i] += 1./(dy[divy]-dy[ind])*prod dpsieta[i] *= fx # print dpsixi[i], dpsieta[i] # Actualiza los valores de las derivadas de la funcion de transformacion dxdxi += x[l,i]*dpsixi[i] dxdeta += x[l,i]*dpsieta[i] dydxi += y[l,i]*dpsixi[i] dydeta += y[l,i]*dpsieta[i] # print dxdxi,dxdeta,dydxi,dydeta for i in range(npe): # Calcula el gradiente para cada funcion de interpolacion grpsi[i,:] = lin.solve(np.array([[dxdxi,dydxi],[dxdeta,dydeta]]), np.array([dpsixi[i],dpsieta[i]])) # Calcula el factor de correccion hs = np.abs(dxdxi*dydeta-dxdeta*dydxi) # Construye la matriz de difusion elemental for i in range(npe): for j in range(npe): elm_dm[i,j] += np.dot(grpsi[i,:],grpsi[j,:])*(hs*wx[nq[0]]*wy[nq[1]]) # Construye los vectores D vecD[i,:] += grpsi[i,:]*(hs*wx[nq[0]]*wy[nq[1]]) # Actualiza la matriz de difusion global Mat_dif[np.meshgrid(C[l,:],C[l,:],indexing='ij')] += mu*elm_dm Mat_dif[np.meshgrid(ng+C[l,:],ng+C[l,:],indexing='ij')] += mu*elm_dm # Extensiones Mat_dif[C[l,:],2*ng+l] -= vecD[:,0] #incorpora Dx en el bloque B Mat_dif[ng+C[l,:],2*ng+l] -= vecD[:,1] #incorpora Dy en el bloque B Mat_dif[2*ng+l,C[l,:]] -= vecD[:,0] #incorpora Dx en el bloque C Mat_dif[2*ng+l,ng+C[l,:]] -= vecD[:,1] #incorpora Dy en el bloque C # Construye el vector del lado derecho vec_b = np.zeros(2*ng+ne) # Implementa las condiciones de frontera for ind in range(ng): if(gfl[ind,0]==1): # Ajusta las entradas al lado derecho vec_b -= Mat_dif[:,ind]*gfl[ind,1] + Mat_dif[:,ng+ind]*gfl[ind,2] # Ajusta el valor de la entrada en b vec_b[ind] = gfl[ind,1] # Coor. x vec_b[ng+ind] = gfl[ind,2] # Coor. y # Ajusta la matriz del problema Mat_dif[:,ind] = np.zeros(2*ng+ne) Mat_dif[:,ng+ind] = np.zeros(2*ng+ne) Mat_dif[ind,:] = np.zeros(2*ng+ne) Mat_dif[ng+ind,:] = np.zeros(2*ng+ne) Mat_dif[ind,ind] = 1. Mat_dif[ng+ind,ng+ind] = 1. return Mat_dif, vec_b

In [6]:
(Mat_dif, vec_b) = matDiff(Ne,Np,ne,npe,ng,x,y,C,gfl,mu)

plt.spy(Mat_dif)
plt.show()

Solución del sistema resultante

Teniendo ya construidos los componentes del sistema a resolver, basta aplicar un esquema de solución razonable para el tipo de matriz obtenida. En este caso se utiliza el solucionador por defecto del módulo linalg. Dado que la matriz construida inicialmente es singular, se omiten los valores de la primera fila y columna - que corresponden a un nodo de frontera y ya se conoce la velocidad en este punto. Esto permite que la matriz obtenida sea invertible.


In [7]:
start = timeit.default_timer()
sols = lin.solve(Mat_dif[1:,1:],vec_b[1:])
print timeit.default_timer()-start

# Recupera las velocidades
solx = np.concatenate([[gfl[0,1]], sols[0:ng-1]])
soly = sols[ng-1:2*ng-1]

# Calcula la magnitud de la velocidad en cada nodo
solm = np.abs(solx**2+soly**2)

# Recupera las presiones
solp = sols[2*ng-1:]


0.240211258681

Gráfica de la solución obtenida

Una vez resuelto el sistema, es posible visualizar el resultado obtenido por medio de una gráfica. En este caso se conoce que la solución analítica corresponde a la linea recta entre las dos condiciones extremas, y se incluye también su gráfica.


In [8]:
# Grafica las fronteras de los elementos
for l in range(ne):
    plt.plot([x[l,0],x[l,Np[0]],x[l,-1],x[l,-Np[0]-1],x[l,0]],
             [y[l,0],y[l,Np[0]],y[l,-1],y[l,-Np[0]-1],y[l,0]],
             'r--')
    plt.hold(True)
    
# Grafica la magnitud del campo de velocidades
for l in range(ne):
    Px = (coords[C[l,:],0]).reshape(Np[1]+1,npe/(Np[1]+1))
    Py = (coords[C[l,:],1]).reshape(Np[1]+1,npe/(Np[1]+1))
    SM = (solm[C[l,:]]).reshape(Np[1]+1,npe/(Np[1]+1))
    plt.pcolor(Px,Py,SM,vmin=min(solm), vmax=max(solm))
    
# Grafica el campo de velocidades
plt.quiver(coords[:,0],coords[:,1],solx,soly,color='w')
plt.title('Campo de Velocidades para el Flujo en la Cavidad')
plt.xlabel('Posicion (coord x)')
plt.ylabel('Posicion (coord y)')
plt.xlim([-1., 1.])
plt.ylim([-1., 1.])
plt.show()

In [9]:
# Grafica el campo de presiones
xpres = np.zeros(ne)
ypres = np.zeros(ne)
for l in range(ne):
    xpres[l] = (x[l,0]+x[l,Np[1]]+x[l,-1]+x[l,-Np[1]-1])/4
    ypres[l] = (y[l,0]+y[l,Np[1]]+y[l,-1]+y[l,-Np[1]-1])/4
    
xpres = xpres.reshape(Ne[1],Ne[0])
ypres = ypres.reshape(Ne[1],Ne[0])
zpres = solp.reshape(Ne[1],Ne[0])

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(xpres,ypres,zpres)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Presion')
#plt.pcolor(xpres,ypres,zpres)
plt.show()

In [ ]: